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SUMMARY 

The aerodynamic design of turtine stator or rotor blades requires 
accurate determination of the velocity distribution on the blades. For 
axial flow turbines a quasi -three -dimensional compressible flow analysis 
has been used successfully for many years. The analysis is based on two 
methods, one to obtain a blade -to -blade variation of velocity, and the 
other to obtain the radial variation in velocity. A combination of the 
methods gives a quasi-three-dimensional compressible flow analysis. Most 
of the calciilations have been incorporated into a program. This 
report gives a description of the quasi -three -dimensional flow analysis 
and the FORTEIAR IV computer program. 


IMEODUCTIOR 

The aerodynamic design of axial flow nozzles and turbine blades re- 
quires accurate determination of the velocity distribution on the blades. 
In the flow analysis, three-dimensional flow effects are of importance. 
Methods for accomplishing this type of analysis for blade designs of 
medium or hi^ solidity have been developed at Lewis Research Center 
(refs. 1 to 4). 

Reference 1 gives a method for determining the variation of velocity 
from hub to tip on a mid-channel stream surface, and reference 2 gives a 
method for determining the blade -to -blade variation velocity. A combina- 
tion of the two methods gives a quasi -three -dimensional flow analysis. 
Weight flow calculations are then based on a two-dimensional integration 
across a passage cross section. Most of the basic calculations needed 
for this analysis have been incorporated into a computer program to facil- 
itate practical turbine design. The program is referred to as CTTD (Com- 
pressor and Turbine Division Turbine Design program) . 

The purpose of this report is to give a complete description of the 
quasi -three -dimensional flow analysis and of the use of the CTTD computer 
program. This technique has been used at Lewis Research Center for 14 
years and the program has long been available to industry for their own 
use in turbine design. The original work of developing the analysis 
method and the CTTD program was done by the staff of the Compressor and 
Turbine Division at RACA Lewis Flight Propulsion Laboratory. 
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SYMBOLS 

a parameter, eq,. (Al) and (A2) 

t parameter, eq. (Al) and (A2) 

C curvature of streamlines on t lade -to -"blade surface, l/ft 
g gravitational constant, ft/sec^ 
m distance along meridional streamline, ft 

n distance along orthogonal to streamline, ft 

nQ distance along orthogonal to streamline on "blade -to -"blade surface 
from suction to pressure surface, ft 

p pressure, Ib/ft 

r radius, ft 

r radius of curvature of streamline in radial-axial plane, ft 

c 

R gas constant, (ft) (l‘b)/(lb) (^R) 

T temperature, deg R 

V absolute velocity, ft/sec 

V weight flow, Ih/sec 

W relative velocity, ft/sec 

z axial coordinate, ft 

a angle of meridional streamline with the axial direction, deg. 

P angle of streamline on hlade to blade surface with the axial 

direction, deg. (see fig. 3) 

X specific heat ratio 

p density, Ib/ft^ 

Od rotational speed, radians/sec 
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Sub scripts 

calc calculated 

cr critical velocity 

e exit 

giv given 

h hull 

i inlet 

m mean 

mid mid- channel 

p pressure 

s suction 

t tip 

d tangential direction 

Superscripts 

' absolute total condition 

” relative total condition 


METHOD AKD ASSUMPTIONS 

The objective of the analysis method is to calculate the quasi-three- 
dimensional velocity distribution satisfying continuity at a given channel 
orthogonal surface (see fig. l). The weight flow may be specified in the 
calculation, or the calciilation may continue until the maxi mum (choking) 
wei^t flow for that channel orthogonal surface is determined- The velocity 
variation on the orthogonal surface is calculated from radial equilibrium 
(ref. l) and stream filament theory (ref. 2). When this is done for several 
orthogonal surfaces, a velocity distribution over the blade surface is deter- 
mined. The velocity distribution can thus be obtained for the guided channel 
formed by the portion of the passage where the blade-to -blade orthogonals 
extend from suction to pressure surface. The guided channel will not cover 
the entire suction surface. To obtain a velocity on the uncovered portion 
of suction surface some method must be used to estimate the location of the 
pertinent stagnation streamline. Since these velocities depend critically 
on the location of the stagnation streamline, unreliable results may be 
obtained on the uncovered portion of the blade. 



-4- 


The "basic simplifying assumptions used in deriving the equations used 

are: 

1. The flow relative to the hlade is steady. 

2 . The fluid is a perfect gas. 

3. The fluid is a nonviscous gas. 

4. The fluid velocity has no radial component. (The projections of 
the streamlines on a radial-axial plane are straight lines parallel to the 
axis. ) 


5. The mid-channel line is a streamline, hereinafter referred to as 
the mid-channel streamline. 

6. (a) The gas has a constant entropy from hub to tip at the mid-channel 
streamline at a fixed axial coordinate. This assumption is used in calculat- 
ing the mid-channel velocity variation (ref. 1). 

(b) The gas has a constant entropy blade -to -blade along an orthogonal 
to the streamlines. This assumption is used in calculating the blade-to- 
blade velocity variation (ref. 2). 

7. A line connecting the midpoints of the hub, mean, and tip orthogonals 
at a fixed axial location in the channel is a straight radial line. This 
assumption is used in calculating the mid-channel velocity variation (ref. 1). 

8. The meridional streamline curvature varies linearly from hub to tip. 
(By assumption 4, meridional streamline curvature would be zero; however, the 
effect of wall curvature on radial equilibrium may be considered. ) 

9. There is free vortex velocity distribution at the inlet to the blade; 
that is, (r"7g) ^ is a constant from hub to tip. This assumption is used in 

calculating the mid-channel velocity variation (ref. 1). 

10. Uniform inlet absolute total temperature, Tj^. This assumption is 
used in calculating the mid-channel velocity variation (ref. 1). 

11. An additional assumption is necessary to calculate the velocity 
variation from blade -to -blade along an orthogonal. An option is provided for 
this assumption in the program. The usual assumption is that there is a linear 
variation of streamline curvature along an orthogonal. An alternate assumption 
is that there is a linear variation of radius of curvature along an orthogonal. 
There is no particular reason why one assumption is preferred over the other. 
For high solidity blading, it makes very little difference which assumption is 
chosen. 

Another possibility is provided in the program. In this case the blade 
surface velocities are calculated as above based on the assumption of linear 
variation of either curvature or radius of curvature. Then for the weight 



-5- 


flow calculation, density and velocity along the orthogonal are computed by 
assuming a linear variation of static pressure along the orthogonal. This 
is not consistent with the assumption of linear variation of either curvature 
or radius of curvature and should be used with caution. 

FLOW CHAUHEL LAYOUT AND ENGINEERIHG DATA 

The initial steps in the design, those involving the development of the 
inlet and outlet velocity diagrams, will not be described (see ref. 5 for 
fiorther information) . It is assumed that these velocity diagrams have been 
obtained, as well as the basic operating conditions of design weight flow, 
number of blades, gas to be used, operating speed and inlet and outlet stagna- 
tion temperature and pressure. From this information an initial blade shape 
can be drawn (ref. 6) . The following steps are followed to obtain a blade 
surface velocity distribution. 

A. Along a section midway between the hub and tip, a cylindrical develop- 
ment of the blade channel should be accurately laid out to scale several times 
actual size. Figure 2 shows a typical blade channel - mean section. (This sec- 
tion need not be at constant radius, but the meridional streamline angle, a, 
should be small enough that cos a is approximately 1.) The mid-channel stream- 
line is then drawn midway between the suction and pressure surfaces, as shown in 
figure 2. This procedure is repeated for the hub and tip of the blade. The 
axial coordinate locations of the hub, mean, and tip sections relative to 
each other must be established. For this, an axial coordinate reference line 
is specified on each of the three blade sections. The relative angular location, 
0, of the hub, mean, and tip sections is not considered; the effect is usually 
neglible . 


B. Any number of axial locations can now be chosen on the three mid- 
channel streamlines. The distance along the mid-channel streamline between 
each axial location shoifLd be measuired for use in specifying loss distribution. 

At any given axial location a curve is drawn through the mid-channel streamline 
from blade -to -blade so as to be orthogonal to each blade surface and to the mid- 
channel streamline, as illustrated in figure 3. This is done at hub, mean, and 
tip. The corresponding orthogonals at the three radial blade stations are posi- 
tioned so that the intersections of the mid-channel streamlines and the ortho- 
gonals are located at the same axial coordinate. In general, the intersections 
of the orthogonals with the suction and pressure sior faces at the three radial 
sections will not be at the same axial coordinates. These three orthogonals 
determine a section through the blade passage over which the velocity variation 
will be determined and across which weight flow will be calculated. The total 
length of the orthogonal, nQ, is required as input for the program. Also the 
blade surface curvatures, Cg and Cp, must be measured at the ends of the or- 
thogonals. These curvatures should be measured very carefully. The angle p 
is the angle that the mid-channel streamline makes with the axis, and is con- 
sidered positive in the direction of rotation, i.e., p is positive if the 
tangential con^jonent of the velocity is in the direction of rotation. 

If the inner and outer walls are not straight, the inner and outer wall 
curvature, l/r^,, should be measured at each station. This curvature is considered 
positive if the wall is concave upwards (fig. 4). This curvature is assumed to vary 
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linearly between hub and tip. The hub and tip radii complete the geometrical 
data required at each station. 

C. Operating conditions must be specified. These include the gas 
specific heat ratio, x, operating speed, OD, and design weight flow per channel, 
w. Also required are and p" at the hub, mean, and tip for each ortho- 

gonal. The inlet relative critical velocity, W^^., can be calculated using 
Vg from the inlet velocity diagram together with the inlet absolute stagna- 
ti6n temperature Tj^. 


W 

cr 



T 


tp** 


( 1 ) 


where 



2or V - (cor)^ 

2rRg 


(r - 1) 


(2) 


By assumptions 9 and 10, T^ and ^ independent of radius, and 

equations (l) and (2) can be used to determine at any point in the 

passage as a function of radius alone. 


The inlet relative total stagnation pressure can be calculated from 


Pi = 


Pi 




r/(r-i) 


(3) 


The exit relative total pressure, p^, can be found in the same manner utiliz- 
ing the exit velocity diagrams. To allow for losses, the difference, p^ - pg, 
is distributed along the length of the mid-channel streamline. Usually a 
linear variation is used. Finally, p" is given by 


P 


tt 



(4) 


D. The information determined above in steps B and C is the information 
for items to and on the Input Data Sheet (fig. 5). One sheet 

is required for each axial station. The calculation at any one station is 
independent of any other. It will be noted that the input sheet is designed 
so that only the numerical values which change from sheet to sheet need to be 
supplied. 

The output from the program includes the blade surface velocities at each 
orthogonal and the corresponding weight flow. These velocity calculations are 
made for the initial estimated value of W/Wg^. at the mid-channel streamline 
for the mean blade section. It is possible for the program to determine the 
velocities corresponding to any desired weight flow (less than choking) or to 
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determine the choking weight flow. If the velocity distribution corresponding 
to the design weight flow or some other operating weight flow is not satisfac- 
tory, the blade design should be altered, and steps A to C repeated to determine 
the new velocity distribution 

DESCRIPTION OF PROGRAM INPUT AND OUTPUT FOR SAMPLE PROBLEM 

The Input Data Sheet is shown in figure 5. The quantities filled in are 
for a san^jle problem. Hie output for this sample problem will be presented 
subs equently . 


General Instructions for Pilling Out Input Data Sheets 

One page must be filled out for each z (axial) station. The basic calcula- 
tion at each station is from a value of (W/W ) (hereinafter called x) to 

mid,m 

the corresponding weight flow velocity distributions. 

These are for problem identification and will be printed as the heading of 
all z stations. Thirty-nine characters are allowed on each line, and they 
a may be alphabetic, numeric, special symbols, or blanks. It is customary to 

n fill out these two lines on the first page only and write OMIT on these two 

d lines on all successive pages. However, if a new heading is desired, a page 
^ which is blank except for the yglue "9" entered in = KR4 must precede the 
(2) page on which the new (l) and ( 4 ) occur. 


(5) JX = if linear variation in 
^ \2/ Iradi 


cmrvature 

radius of curvature/ 


is to be assumed in 


the calculation of the velocity distributions across the orthogonals. 


® JY = Jxj 


linear variation in 


the calculation of the wei^t flow. 


/same 

\stati 


same as JX 
static pressure 


is to be used in 


® JZ = 1 if only the calculation for the one value of = x is wanted 


In this case. 


= w . 
giv 


is not used. 


JZ = 2 if choke conditions are wanted. (Solve for the value of x which 


gives maximum value of w , .) In this case w . is not used. 

calc giv 


If there is any value of x for which = '^eiv equal to 

there are always two such values. ^ 


JZ = 3 if the lesser of these two solutions is wanted, (subsonic) 

JZ = 4 if the greater of these two solutions is wanted, (supersonic) 
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3 

14) 


X, ratio of specific heats 
CO, wheel speed, (radians per second) 
w weight flow, PER BLADE, (^/sec) 
First page. 

Fill out (l^ (l4) and (^) = KR4 = 1. 


On any successive page, if values to Tdc entered in 

ferideSlcal to) in 0 - 0 of the immediately preceding 

J) 


V 

page, KR4 = 


If KR4 = 0, write OMIT in (12 


= KR4 = 9 may l^entered on an otherwise blank page. See explana- 
tion under and 


(‘" 1 . 

W. 


reciprocal of radius of curvature at hub, (ft“^) 
reciprocal of radius of curvature at tip, (ft“^) 


XI) 


radius at hub, (ft) 
r^, radius at tip, (ft) 

First page. 

Fill out 0 - 0 and (T) = KR5 = 1. 




On any successive page, if values to be entered in (18 

( differ from 
are identical tq, 

page, KR5 = 

, (0 , (^ Uq, length of the orthogonal between suction and pressure 
surfaces at hub, mean, tip, (ft) 


those in (^ (1^ on the immediately preceding 

If KR5 = 0, write OMIT in 0 0 • 


20) , (2^ , (^ Cg, curvature at intersection of orthogonal with suction 

surface at hub, mean, tip, (ft“^) 


, ^26) , (31) C , curvature at intersection of orthogonal with pressure 
surface at hub, mean, tip, (ft“^) 


both must 
be posi- 
tive (non- 
zero) if 
JX = 2 


(B) 


, , (g) relative critical velocity at hub, mean, tip, (ft/sec) 

, P"’'^cr 'Weight flow parameter at hub, mean, tip (^/ft^sec) 
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First page. 


Fill in @ • 
If values for 
write OMIT in 
them and 




and -* repeat those of -* , 

and = KR6 =1. If values differ, enter 


On any successive page, if all values to he entered in -* are 

identical to th^e in -♦ of the immediately preceding page, 

write OMIT in -♦ (3^ and KR6 = 0. If not, follow instructions 

given above for the first page. 


t P, relative (rotor) or absolute (stator) flow angle measured from the 
axis, (degrees), positive in the direction of rotation, mid-channel 
only. 


First page. 





ijB constMt from hub to tip, 
@ and @ = KR7 = 1. 


fill in 


only, write OMIT in 


If P is specified Qftly at hub, mea^and tip, fill in 
, write OMIT in and ^ = KR7 = 3. 


and 


If P y^riation^rom hub to tip is specified (nine values), fill in 
0) - @ and @ = KR7 = 9. 


On any successive page, if all values in ^54} ^^ (42^ are id^xtical to 
those of the preceding page, write OMIT in (3^ and = KR7 

=0. If not, follow instructions given above for first page. 


An ei^t digit numerical code to identify this page uniquely is entered 
here and will be printed out preceding the output for this z station. 

(Can be used to code for engineer, case and z-station number.) 

I , critical velocity ratio, (referred to as x immediately below) . 
niid,m 



Computation of ''^ga2c the velocity distributions) for this value 

of X will be performed first. If; 

© = JZ = 1, calculation stops, and input data for next station is read. 

© = JZ = 2, X is increased automatically and calculations are continued 
until choking weight flow has been found. 

© = JZ = 3, X is modified until subsonic solution is found. (Point at 
which 'W'calc ~ ^giv^’ ®^®ss as close as possible. 

© = JZ = 4, X is modified until supersonic solution is found. Guess 
LOW, 
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Description of Computer Output 


An example of the output obtained is given in table I. The first output 
is a listing of the information on the Input Data Sheet. In the sample output 
these are identified with the numbers corresponding to the numbers on the Input 
Data Sheet. If nine values of P (equally spaced from hub to tip) are not sup- 
plied, the missing values are calculated and printed with the input values. The 
remaining output is the calculated quantities. The first line is and 


Wp/Wmid^ each at hub, mean, and tip. Next are the values of the parameters a 

and b of equation (A2), which are labelled LITTLE A and LITTLE B. IK = 1 is 
the hub and IK = 9 is the tip. ALPHA is equal to a(r.j. - rj^)/l6. M(IK) is 

the numerical approximation to j a(^)d^ and R(IK) is be PH^ PM, 

p I 




and FT are the numerical approximations to 


b(|)e 


d^ at the 


m 


hub, mean, and tip, respectively. QH, QM, and QT are the numerical approxi- 


r 


a(0d? 


mations to e . The relative velocities ratioed to the relative crit- 

ical velocity at the hub, mean, and tip, for suction surface, mid-channel and 
pressure surface are followed by the values of A(I,K) ( pW/p"W(,j.) . The first 


line is at the hub at 6 equal intervals from the suction surface to the pressure 
surface. The second and third line of A(IK) are at the mean and tip, respec- 
tively. On the next line R is the number of the iteration for this orthogonal, 
X is the value of (wA^cr)mid m this iteration, WT FIOW CALC is the 

corresponding calculated weight flow, and D is the measure of error 

^ ■ ^calc/^giv' more than one iteration is done, the values down to PH, PM, 

FT, etc. , do not change, hence for the second and following iterations only the 
lines following this are printed. A maximum of 5 iterations is performed. If 
the convergence criterion is not satisfied, the statement "NLIMIT (5) HAS BEER 
REACHED" will be printed. 


The message CALCULATION OF A (IK) AT ST. NO. 152+5 IN CTTD MP ASKS FOR 
LOG OF REG. REMAINING CALCULATIONS THIS ITERATION THEREFORE INVALID. is usually 
caused by an error in the input cards. If this is not the case, the blade design 
must be changed. 


CTTD COMPUTER PROGRAM 

The program consists of the main program and two subroutines, PABC and 
VSUBX. PABC calculates the coefficients A, B, and C of the parabola 
y = Ax^ + Bx + C passing through three given points. VSUBX calculates the 
ratio of the velocity at any point on an orthogonal to the mid-channel velocity, 
using equation (A4) or (A5). 
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C CTTO MAIN PROGRAM 

DIMENSION ENO< 3),CS(3) ,CPm tWCRC3) , PHCR(3 ) t VS ( 3 I , VP( 3 ) , WM ( 3 ) , 

1S( 3)» BETA(9},ALFt 9|,EM(9) » C APN ( 91 , A( 9) » B( 9 ) » MS { 3 ) t WP( 3 ) 

001 EQUIVALENCEIWSmtWSHI t(WS(2) fWSM), ( WS ( 3 ) » WST) , ( WP{ 1 ) , WPH) , ( WP( 2 
!)• WPM),(WP(3),WPT) 

DATA NL IM,TLIM/5,.00l0/ 

400 FORMAT ( 40HI ) 

401 FORMAT I 40H ) 

402 FORMAT (Till 

404 FORMAT ( 3E9.5I 

406 FORMAT (4E9.5) 

408 FORMAT (5E9.5) 

410 FORMAT ( IHO, 2 1 4, E9. 5 » 

412 FORMAT (llHO GAMMA = ,E12.5,11H OMEGA = ,E12.5,12H WGIVEN = ♦ 
1E12.5) 

414 FORMAT ( I06H0 N ZERO C SUB S C SUB P 

1 W SUB CR RHO W CR I /R SUB C R) 

416 FORMAT ( lOHO HUB t7El5.5| 

417 FORMAT < lOH MEAN »5E15.5) 

410 FORMAT ( lOH TIP t7El5.5) 

420 FORMAT ( 1 7H0 BETAI 

422 FORMAT ( IHO, I 7, 1 4, E 1 5. 5, 3 14 ,I13,3Il///» 

430 F0RMATI/6H0 VSH=,E 12. 5, 6H VSM= E12.5,6H VST= , E12 . 5 t6H VPH=,E12. 
15*6H VPM=,E12.5,6H VPT=, E12.5 ,/45H0 IK LITTLE A LITTLE 

2B ALPHA /) 

432 FORMAT ( I4,3F15.5) 

434 FORMAT ( 31H0 IK MIIK) NIIKJ /) 

436 FORMAT ( 14,2E15.5I 

438 FORMAT ( 6H0 PH=,E12.5,6H PM= E12.5»6H PT=,E12.5,6H QH=,E12. 

15, 6H 0M = ,E12.5,6H 0T=,E12.5) 

440 FORMAT C//9H0 WSH = ,E12.5,9H WMH = ,E12.5,9H WPH = ,E12.5,/9 
IH WSM = ,E12.5,9H WMM = ,E12.5,9H WPM = ,E12.5,/9H WST = , 
2E12.5,9H WMT = ,E12.5,9H WPT = ,E12.5//) 

442 FORMAT ( 8H A UK ) =, 9E 12. 5 ) 

444 FORMAT < 7H0 N = ,I1,16H X * WMM/HCR = ,E12.5,18H WT FLOW CALC 
1= ,E12.5,7H D = ,E12.5) 

006 READ! 5,400) 

008 READ (5,401) 

010 READ (5,402)JX, JY,JZ*KR4,KR5,KR6,KR7 

012 IF(KR4) 25,14,25 

014 IFIKR5) 30,16,30 

016 IF(KR6) 18,20,18 

018 IF(KR6-1) 42,34,42 

020 IF(KR7) 22,66,22 

022 IF(KR7-1) 24,48,24 

024 IF(KR7-3) 64,56,64 

025 IF(KR4-9)26,6,26 

026 READ (5,404)GAMMA,0MEGA,WGI V 
028 GO TO 14 

030 REAO (5,406)RCH,RCT,RH,RT 
032 GO TO 16 
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034 READ { 5, 408JEN0( l),CS( 1) ,CP( 1) tWCR ( U f PWCR( U 
036 00 38 K=lt2 

EN0{K+l)=EN0(K) 

CS (K>l)=CS (K) 

CP(K+1I=CP(K) 

WCR(K»l)=WCR(K) 

038 PWCR(K+l)=PWCR(KI 
040 GO TO 20 
042 DO 44 1 = 1,3 

044 READ (5,408)EN0( n,CS(n,CPII),WCR(I),PWCRjn 
046 GO TO 20 

048 READ (5,404JBETA< U 

050 00 52 IK = l,8 

052 8ETA( IK4-1)=BETA( IK) 

054 GO TO 66 

56 READ (5,404)BETA( n ,BETA(5),BETA(9) 

5 8 DO 60 K=l,5,4 

TEMP=(BETAIK+4)-BETA(K) ) /4.0 

59 DO 60 J= 1,3 
KJ=K+J 

BETA(KJ)=BETA(KJ-l)+TEMP 

60 CONTINUE 
062 GO TO 66 

64 READ (5,404) (BETA(IK),IK=1,9) 

066 READ (5,410) 101, 102, X 
068 WRITE (6,400) 

070 WRITE (6,401) 

072 WRITE (6,412)GAMMA,0ME6A,WGIV 
074 WRITE (6,414) 

076 WRITE (6,416)EN0( l),CS(l),CP(l),WCR(l),PWCR(l) ,RCH,RH 
078 WRITE (6,417)EN0(2),CS( 2),CP(2) ,WCR(2) ,PWCR(2) 

080 WRITE (6,418)EN0( 3),CS( 3) ,CP(3) ,WCR(3) ,PWCR(3) ,RCT,RT 
082 WRITE (6,420) 

084 WRITE (6,416)BETA( 1),BETA(2),BETA(3) 

086 WRITE (6,417)BETA(4),BETA(5),BETA(6) 

088 WRITE (6,418)BETA(7),BETA(8),BETA(9) 

090 WRITE (6,422) 101,102, X,JX,JY,JZ,KR4,KR5,KR6,KR7 
100 PM = 0.0 
QM = 1.0 
MPSET = 1 
NC =1 
KSTAR = 1 
N =0 

XOLO = 0.0 
TWOOM = 2.040MEGA 
Cl = GAMMA- 1.0 
C2 = GAMMA/Cl 
C3 = C1/(GAMMA + 1.0) 

C4 = 1.0/GAMMA 
C5 = 1.0/Cl 
OR » RT-RH 
0R6 * OR/6.0 

0R16 = DR/16.0 

DR24 = OR/24. 0 
DRC = RCT-RCH 
P=0.0 



- 13 - 


lOl DO 104 1 = 1,3 

CALL VSUBX (ENOU)tCSm,CPII),P,JX,VX) 

104 vsc n = vx 

106 P = 1.0 

108 00 110 I = 1,3 

CALL VSUBX CENOCII,CS(n,CPni,P,JX,VX) 

110 VPC II = vx 

112 WRITE C6,430IVSm,VSf 2I,VSC3I,VPU),VPI2I,VPI3) 

114 P = 0.0 
116 00 120 IK=1,9 
RK =RH-»-P»DR 

RCK =RCH*P*1»C 
SNB =SINIBETAI IKI*. 01745329) 

B(IK) =TWOON*SNB 

SNB2 =SNB*SNB 

CSB2 =l.0-SNB2 

AK =(CSB2*RCK)-(SNB2/RKI 

ALFnKl=AK*0R16 

118 WRITE (6,432)IK,AK,BUK),ALF(IK) 

120 P = P+.125 

122 EM(4l = -ALF(5)-ALFI4I 

EMI 3) = EMI4I-ALF(4)-ALF(3) 

EMI 21 = EM(3)-ALF< 3I-ALFI2I 
EMU) = EM(2)-ALF(2)-ALFI 1) 

EMI6) = ♦ALF(5)+ALF(6) 

EMIT) = EMI6)*ALF(6I+ALFI7) 

EMI8) * EMI7)+ALF( 7M-ALFI8) 

EMI9) = EM(8)*ALF( 8H-ALFI9) 

EMI 5) =0.0 

123 WRITE 16,434) 

124 DO 126 IK =1,9 

CAPNl IK)=BIIK)/EXPIEMIIK|) 

126 WRITE I6,436)IK,EMIIK),CAPNIIK) 

128 PH=-0R24»ICAPMI1)*CAPNI5I+2.0*CAPNI3)+4.0*ICAPNI2I+CAPNI4) ) ) 
PT=^DR24»ICAPNI5H'CAPNI9)<-2.0*CAPNI7)^4.0*ICAPNI6)+CAPNI8) ) ) 
QH =EXPIEMI1)I 
OT =EXPIEMI9)I 

130 WRITE I6,438)PH,PM,PT,QH,QM,QT 
140 WMIl) =QH*IX*WCRI2)-PH)/WCRI1I 
WMI3) =QT*IX*WCRI2)-PT)/WCRI3) 

UM f 7 1 = Y 

WSH =WMI1)*VSII) 

WST =WMI3I«VSI3) 

WPH =WMI 1)*VPI 1) 

WPT =WMI3)*VPI3) 

WSM =X*VSI2) 

WPM =X*VPI2) 

142 WRITE I6,440)WSH,WMI1I ,WPH,WSM,WMI2) ,WPM,WST,WM13) ,WPT 
144 DO 166 1 = 1,3 
146 SIGA =0.0 
P=0.0 

ENOC =EN0I I I 

148 GO TO I 150, 150, 158), JY 
150 CPC =CPII) 

CSC =CSIII 
WMC =WMII) 
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152 DO 154 IK = 1,9 

CALL VSU0X(ENOCtCSC,CPC,P, JX,VX) 

IK =VX*WMC 
TEMPO * l.-C3*ZK*ZK 
IF(TEMP0»1153, 1153,153 
1153 TEMPO * ABSf TEMPO) 

WRITE(6, 1234) 

1234 F0RMAT(56H CALCULATION OF A(IK) AT ST. NO. 152+5 IN CTTO MP ASKS 
12341, 60HFOR LOG OF NEG. REMAINING CALCULATIONS THIS ITERATION THERE , 
12342 13HFORE INVALID. ) 

153 A(IK) = ( TEHPO**C5)*ZK 
SIGA s$lGA+A(IK) 

154 P =P+.125 
156 GO TO 164 

158 ELS =( l.0-C3*WS( I)*WS( I ) )**C2 
ELP =( l.O-C3*WP( I )*WP( I ) )**C2 
160 00162 IK»l,9 

CAPJ =ELS+P»(ELP-ELS) 

TEMP =CAPJ**C4 

A( IK) =TEMP*SQRT((l.0-CAPJ/TEMP)/C3) 

SIGA sSIGA+AI IK) 

162 P =P+.125 

164 WRITE (6,442MA( IK),IK=l,9) 

TMP *.5*(A( 1)+A(9) ) 

166 SID =( SIGA-TMP)*(ENOC/8.0) 

168 N »N + l 

W *0R64| Sm*PWCRm+4.0*SI2)*PWCR(2) + S{3)*PWCR(3) ) 

0 =1.0-{ W/WGI V) 

170 WRITE (6,444)N,X,W,0 

172 GO TO ( 10, 174, 174, 174), JZ 
174 IF (NLIH-N) 176,176,180 
176 WRITE I6,500)NLIM 

500 FORMAT ( 13H0 NLIMIT (,I3,18H) HAS BEEN REACHED ) 

178 GO TO 10 

180 GO TO (182,186, 198, 198),JZ 
182 WRITE (6,501) 

501 FORMAT ( 19H ERROR STOP AT 180 ) 

184 GO TO 10 

186 GO TO (188,192,256),KSTAR 

188 KSTAR =2 

190 GO TO 210 

192 KSTAR =3 

194 GO TO 240 

198 IFl ABS{0I“TLIM) 10,200,200 
200 GO TO (202, 230, 242, 256), NC 
202 NC=2 

204 IF( ABS(D)-. 05)206, 210, 210 
206 DX =.05 
208 GO TO 212 
210 DX =.20 
212 X1=X 
WI=W 

01 = 0 

214 IF(D.GT.O.O) X=X+OX 
216 IF(D.LT.O.O) X=X-DX 
226 GO TO 140 
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230 NC=3 

232 IF{0*DU238f 234*240 
234 WRITE (6*503) 

503 FORMAT ( 19H ERROR STOP AT 232 I 
236 GO TO 10 

238 OX =.540X 
MPSET *2 

240 X2=X 
W2=W 
02=0 

241 GO TO 214 

242 IF( 0*021 248, 244* 250 
244 WRITE (6,504) 

504 FORMAT C 19H ERROR STOP AT 242 ) 

246 GO TO 10 

248 MPSET =2 

250 GO TO (252,256), MPSET 
252 X1=X2 
W1=W2 
01=02 
X2=X 
W2=W 
02=0 

254 GO TO 214 
256 X3=X 
W3=W 
03=0 

CALL PABC(X1,X2,X3,W1,W2,W3*APAB,BPAB,CPAB) 

260 GO TO (262,266,290,290) *J2 ' 

262 WRITE (6,505) 

505 FORMAT ( 19H ERROR STOP AT 260 ) 

264 GO TO 10 

266 X =-BPAB/(APA8+APA8) I 

268 IF( .OOl-ABS(X-XOLO) 1270,270,10 

270 XOLO=X I 

272 IF(W1-W2) 274,280,280 j 

274 IF(W1-W3)276,140,140 
276 X1=X3 
W1=W3 
01=03 

278 GO TO 140 

280 IF(W2-W3)282,140, 140 

282 X2=X3 i 

W2=W3 
02=03 

284 GO TO 140 

290 OISC=BPAB**2-4.0*APAB*(CPAB-WGI V) 

292 IF (DISC) 294,298,298 
294 WRITE (6,506) 

506 F0RMAT(55H PABC FIT GIVES NEGATIVE DISCRIMINANT. PROBABLY CHOKED.) i 
296 GO TO 10 

298 XPL=(-BPAB+SQRT(OISC))/(APAB»APAB) i 

XMI=I-BPAB-SQRT(OISC) )/(APAB+APAB) 

300 GO TO (302, 302, 306*308), JZ j 

302 WRITE (6,507) 

507 FORMAT ( 19H ERROR STOP AT 300 ) I 
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304 GO TO 10 

306 X = AMINUXPL,XMI ) 

GO TO 316 

308 X = AMAXl(XPL,XMI ) 

316 IF ( ABS{0U-ABS(D2) )318,320,320 
318 IF I ABSI D2)-ABStD3M328,322,322 
320 IF ( ABSl Dl)-ABS(D3n328,326,326 
322 X2=X3 
W2=W3 
02=03 

324 GO TO 328 
326 Xl=X3 
W1=W3 
01 = 03 
328 NC=4 
330 GO TO 140 
END 


C SUBROUTINE PABC 

SUBROUTINE PABC { X I , X2»X3, Wl , W2,H3, A, B ,C » 

C9 = XI+X2 

Cll=Xl*Xl 

C15=X3-XI 

C18= (M2-WII/(X2-X1) 

A= { C 1 5*C 1 8-W3+W 1 ) / ( C I 54C 9-X3*X3+C 1 1 ) 

B=C18-C9»A 

C=Wl-Xl*B-Cll*A 

RETURN 

END 


C SUBROUTINE VSUBX 

SUBROUTINE VSUB X I EN.CS.CP.P. JX.VX) 

700 GO TO ( 7C2,706), JX 

702 VX=EXP(EN*ICS*I .5-P ) +.5* (CS-CP » ♦ ( P*P-.25 )) ) 

704 GO TO 714 

706 IF ICS-CPI 712.708,712 
708 VX = EXP( EN*CS*( .5-P n 
710 GO TO 714 

712 VX= (2.0*( IP*(CS-CP)+CP)/(CS+CPI n**(EN4CS*CP/ICP-CS) ) 
714 RETURN 
END 



- 
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List of Program Variables in Main Program 


A 

array for nine values of ■) 

P Wqj. 

- 

AK 

a 


ALF 

(rt - rh) 

array for nine values of a — 

16 


APAB 

coefficient of calculated by subroutine PABC 


B 

array for nine values of b 


BETA 

array, input, p 


BPAB 

coefficient of x calculated by subroutine PABC 


Cl 

r - 1 


C2 

r 

r - 1 



r - 1 



r + i 


C 4 

1 

r 


C5 

1 H 

H 


CAPJ 



CAPE 

, -M 
array, be 


CP 

array, input, Cp 


CPAB 

constant coefficient calculated by PABC 


CPC 

ten 5 )orary storage for current value of CP 


CS 

array, input, Cg 

- 

CSB2 

cos^P 


CSC 

ten 5 )orary storage for current value of CS 


D 

^ Wcaic 
V . 

giv 
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D1,D2,D3 

DISC 

DR 

DR16 

DR24 

DR6 

DRC 

DX 


ELP 

ELS 

EM 

ENO 

ENOC 

GAMMA 

I 

IDl 

ID2 

IK 

J 

JX 

JY 

JZ 


ten^orary storage for various values of D 
temporary storage 


^t - 


At - 
At - 

(l/r,)^ - (Vr.)^ 

arbitrary change in 
of weight flow 

P" 

Ps 


(■ 


X lvalue of 



for next calculation 


J '^r 

a dC (Trapezoidal rifle is used) 
array, input, n^ 

temporary storage for current value of ENO 
input, T 

index of most loops executed 3 times (for hub, mean, tip) 

input, z station code, 4-digit number supplied as input and 
printed out for identification 

input, z station code, 4 -digit number supplied as input and 
printed out for identification 

index of most loops executed 9 times 

index of an input storage loop 

input 

input 

input 
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K index for 2 loops controlling storage of input 

KJ ten 5 )orary storage K + J 


KR4,KR5,KR6,KR7 inpiat, read control switches 


KSTAE 


MEBET 


N 


NC 

NLIM 


OMEGA 

P 

m 

PM 

PT 

PWCR 

QH 

QM 

QT 

RCH 

RCK 


switch-controls branching for successive values of x, 
if choke solution is wanted. Initialized by program 
and automatically stepped. 

switch-controls branch to continue calculations or call 
subroutine PARC. Initialized by program and altered 
as a result of calculations. 

counter for number of times weight flow is calculated. 
When N = NLIM calculation for that page of input 
is stopped and a message written at end of output. 

switch-controls branching for successive values of x 
when subsonic or supersonic solution is wanted 

constant (supplied in DATA statement) limiting number of 
times weight flow calcialation can occur for 1 sheet of 
input data 

CD 


be~®^ dr (Simpson’s rule is used) 


be”^ dr (SiB5>son’s rule is used) 

input, p"Wgj. 
e“'" at hub 
1.0 

e®^ at tip 
input, 

current value of l/r^. 


input, 

n 

“0 



array. 
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RCT 

input, (l/r ) 

^ t 


RH 

Input, 


RK 

current value of r 

- 

RT 

input, r^ 


S 

array, values of { iiR— ■ dn at hub, mean, and tip 

UO ^ cr 


SIGA 

tenporary storage, ^A(IK) 


SHB2 

sin^P 


SUB 

sin p 


TEMP 

temporary storage 


TEMPO 

temporary storage 


TLIM 

constant (supplied in DATA statement). Calculation is terminated 
^calc 

if TLIM > 1 - whenever subsonic or supersonic solution 

^giv 

is sought. 


TMP 

temporary storage 


TWOOM 

20) 


VP 

array, three values of velocity (calculated in subroutine VSUBX) 
for pressure surface at hub, mean, and tip 


VS 

array, three values of velocity (calculated in subroutine VSUEX) 
for suction surface at hub, mean, and tip 


VX 

ratio of velocity to mid-channel velocity 


W 

/ W \ 

weight flow calculated for current value of — I 

mid,m 

• 

W1,W2,W3 

temporary storage for various values of weight flow 


WCR 

array, input, 


WGIV 

input, 


WM 

w 

array, mid-channel values of rj — at hub, mean, and tip 

^cr 




- 21 - 


WMC 

WP 

WPH 

vm 

WPT 

W5 

WSH 

WSM 

WST 

X 


X1,X2,X3 

XMI 

XOLD 

XPL 

ZK 


current value of W! 

W 

array, — on pressure surface at hub, mean, and tip 
^cr 

WP(l) 

WP(2) 

WP(3) 

w 

array ^ tj — on suction surface at hub, mean, and tip 
^cr 


WS(1) 

WS(2) 

WS(3) 


input. 


eJLso the current value of 



for which the 


weight flow is calculated 

ten 5 )orary storage for various values of x 

the lesser of the two values of x at which the line y = 

intersects the parabola y = APAB(x^) + BPAB(x) + CPAB 

temporary storage for a value of x when choking weight flow is 
so;aght 

the larger of the two values of x at which the line y = 
intersects the parabola y = APAB(x^) + BPAB(x) + CPAB 


W 

current value of ^ — at any point 

"cr 
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APPEKDIX - MALYTICAL EQUATIONS 


In reference 1 a differential equation is derived for the variation 
of velocities in a radial direction based on the assun5)tion of a mid- 
channel stream surface (axial symmetry) which is radial at any fixed axial 

position (i.e.,^=0 on the midchannel stream surface). This equation is 

•^ = aW - b I 


where 


r cos a 


V • 2(jd ^ 

b = sin 6 1 + tan a , 

^ Vcos a dm J 


The coordinate system is shown in figure 6. 

With the assumption that a is sufficiently small so that cos a 
is nearly 1 and tan a is nearly zero, we can set dn = dr to obtain 


= aW - b 
dr 


where 


b = 2 od sin p 

An analytical solution of equation (A2) is 


f 

J^m 




- _r 


b(|)e 


-IJ 


for the midchannel stream sheet. 




^mld m computed from the input values 


mid,m 


44) and 


With equation (A3) the midchannel velocities at the hub and tip can then 
be calculated from the value of jji* computing the integration in 
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eq.'oation (A3), the interval to is divided into eight equal 

intervals. Since the integration is started at the mean radius, only- 
four integration steps are required in each direction. The integration 
is done numerically using the trapezoidal rule. In the case where 3 
is given as input only at hub, mean, and tip (rather than at all ninp. 
stations) linear interpolation is used to determine the values which are 

not given. The meridional streamline curvature, s^, is assumed to vary 
linearly from hub to tip. 


The ratios of surface velocity to midchannel velocity. 


can be computed based on an ass-umption of linear -variation of 


either streamline curvature or streamline radius of curvature. These 


equations are derived in the next section. If linear variation of curva- 



where n is the distance along the orthogonal from the blade suction s-ur- 
face. Thus, the velocity can be determined at ar^ point on the orthogonal 
from (A4) . If, on the other hand, linear variation of radius of curvat-ure 


is assumed 


W 

mxd 



(A5) 


By substituting n = 0, and n = in equation (A4) or equation (AS), 
and losing the results from equation (A3)^ the blade surface velocity at 
hub, mean and tip can be confuted. 

With the velocities obtained, the weight flow past the orthogonal 
can be computed. This is calculated from 



The inner integral is confuted at hub, mean, and tip using eight equal 
intervals from suction to pressure surface, and using the trapezoidal rule. 
If JY = 1 or 2 (linear curvature or radius of curvature variation) the 
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W 

critical velocity ratios, , are computed using the input value of W, 


cr 


and equations (A4) or (A5), and (A3). The weight flow parameter, -■ , 

P Wcr 

is then calculated from 

1 


pW _ 

p'V, “ 


'cr 


1 - 


r - 1 / w \ 
r + 1 


r-i 


w 


w 


(A7) 


cr 


If JY = 3 (linear variation of static pressiire), then the blade surface 
static pressures are calciolated based on the svirface velocities 




/ H Y 


JL . 

r-1 


and then pressures along the orthogonal are calciilated by 


then 


L_£s £sl 

P P n^\^P P 


(A8) 


(AS) 


M. 


p'X 


cr 


1 




b'7 ' 

y + 1 
r-1 

;-W'J 

I 




J 


1/2 


(AlO) 


After computing the inner integral in equation (A6) at hub, mean, 
and tip^ the outer integral is approximated by Simpson's Rule. If 
JZ = 1, no further calculations are made. 


If JZ = 2^ further calculations are made with new values of 

(—] 

to get the weight flow, w, as a function of I / . , 

\ cry mid, m 

W \ 


±L OL 

i^] 

\ cry mid, m 


Using a parabolic approximation, the value of 


'crl-: 


weight flow can be estimated. When two successive estimates of I ^ 


giving maximum 


mid^m 


w \ 


^^Aiid^m 

differ by less than 0.001^ the computation is stopped. The maximiun weight 
flow is normally determined within 5 iterations . 
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If JZ = 3 or 4 further calculations are made to determine a 
value of (w/w ) . . which will give a weight flow, w, which is 
close to the input value, w . • If w is less than choking 

giv giv 

weight flow, there are two values of (w/W ) . ^ which will give 

cr mid,m 

w . = w . - If JZ = 3 the smaller, or subsonic, value of 

calc giv ^ ^ 

(w/w^r)jjiid ni found. If JZ = 4 the larger, or supersonic, 

solution will be found. When w . - w , <0.001 w , the calcula- 

I giv calc giv 

tions are stopped. A solution is normally found within five iterations. 


If w . is larger than the choking weight flow no solution exists, 
giv ^ 

In this case calculations are made for five values of (w/w ) . ^ 

' cr mid,m 

and the choking weight flow can usually be estimated from these values. 


Blade to Blade Velocity Variation 

The method of calculating blade surface velocities from mid-channel 
velocities is based on reference 2. The assumptions that the flow is 
steady relative to the blade, nonviscous, and isentropic along the blade 
to blade streamline orthogonal yield 


dn " r^ 


(All) 


This is equation (l) of reference 2 and can be derived from the force 
equation. 

If it is assumed that the curvature, C = l/r^, varies linearly 


along the orthogonal, the equation can be integrated to obtain 
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2(Cp - C3, 


fC + C 


which is equation (7) of reference 2. Since the cvirvature is assumed to 
vary linearly, 

C = C + (C - C )— (Ai 

s ' p s'n^ 

When equation (A13) is substituted in equation (A12), equation (A4) is 
obtained. 

If it is assumed that the radius of curvature varies linearly, then 
“■c ' 

Using this in equation (All) we have 


Integrating equation (A15) from the mid-channel gives 

(r ) + j(r ) - (r ) I 

c s 1 c p ' c s I 

W n .2 

log r~ = (r ) ■- Vr T 7~r~77~^ , X In 

mid cp cs (^c^s r"*c^p - ("'c^s ir 

Substituting ^ into equation (A16) and solving for ^ 


equation 


(AS) is obtained when C 4 C • For the special case where C_ = C„, r 

s p s p c 

is a constant in equation (All) ^ and an integration gives 


mid 


•■c ® % 


which is equivalent to actuation (AS) when Cg = Cp. 
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TABLE I 


% 


Q SAWL€ PROBLEM 
0 NEGATIVE REACTION TURBINE 

^GANMA s 0.14000E 01© OMEGA = 0-48800E 0? 


HUB 

MEAN 

TIP 


©WGI V 


» 0.59170E 00 

W SUB CR 


© 0.42080E-01 i 

g 0. 15710E 

02 i 

0.88200E 

01 1 

j|o,94970E 

03 © 

U0.545OOC-OI { 

5 0.14160E 

02 j 

[I 0.84200E 

01 { 

5 O- 9574 OE 

03 © 

0O.715OOE-O1 ^ 

g 0-12000E 

02 ^ 

5 0.63000E 

01 ( 

g 0.96660E 

03 0 


RHO W CR 1/R ' 

0-47160E 02 © 0 
0.49520E 02 
0.52540E 02 


0 . 


HUB 

@-0-50000E 

01 -0.78250E 

01 

-0.106B0E 

02 

MEAN 

-0.13475E 

02 @-0. 1&300E 

02 

-0-193B0E 

02 

TIP 

-0.22400E 

02 -0.25450E 

02 

©-0.26500E 

02 

1000 

3 0.80000E 00 1 1 

3 

^U33^ 


@ 


© © o 


R 

© 0. 87500E 00 
@ 0. 12500E 01 


VSH* 0.13422E 01 VSM= 0.14145E 01 VST= 0.14594E 01 VPH= 0-B0106E 00 VPM= 0. 76449E 00 VPT= 0. 75858E 00 


IK 

LITTLE A 

LITTLE 8 


ALPHA 

1 

-0,86813E-0Z 

-0.85064E 

02 

-0-20347E-03 

2 

-0.20107E-01 

-0-13288E 

03 

-0.47126E-03 

3 

-0,35256€-0l 

-0. 18037E 

03 

-0.82632E-03 

4 

-0.53463E-01 

-0.22743E 

03 

-0.12530E-02 

5 

-0. 74140 E-01 

*0.27393E 

03 

-0,l7377E-02 

6 

-0.98961E-01 

-0, 32339E 

03 

-0.23194E-02 

7 

-0.12559E 00 

-0.37192E 

03 

-0.29435E-02 

8 

-0.15349E 00 

-0.41941E 

03 

-0.35973E-O2 

9 

-0.18214E 00 

-0.46571E 

03 

*0- 42690E-02 

IK 

M( IK) 

N(IK) 



I 

0.70424E-02 

-0.84467E 

02 


2 

0.63676E* 02 

-0, 13204E 

03 


3 

0.50701E-02 

-0.17946E 

03 


4 

0.29907E-02 

“0.22675E 

03 


5 

0. 

-0.27393E 

03 


6 

-0.40571E-02 

-0,324 TOE 

03 


7 

-0.93200E-02 

-0.37541E 

03 


8 

-0. 15861E-01 

-0.42612E 

03 


9 

-0.237276-01 

-0,47A89E 

03 



PH= 0.33632E 02 PM= 0, 


PT=-0.70389E 02 


0H= 0.10071E 01 0M= 0. lOOOOE 01 


QT^ 0.97655E 00 


HSH = 

0.10422E 

01 

WMH = 

0.77692E 

00 

MPH = 

0-62204E 

00 

HSM ^ 

0- H316E 

01 

VMM = 

o.Booooe 

00 

*<PM = 

0.61159E 

00 

WST = 

0. 12331E 

01 

WMT = 

0.84492E 

00 

WPT = 

0.64103E 

00 


AI IK) = 

0.63259E 

00 

0.63282E 

00 

0.5 2495E 

00 

0-61191E 

00 

0.59588E 

00 

O,57041E 

00 

0-8608 7E 

00 

0-543126 

00 

0.52655E 

00 

AtlK>= 

0.62095E 

00 

0.b3325E 

00 

0.63139E 

00 

0.62024E 

00 

0.60343E 

00 

0.SB3S5E 

00 

0.56239E 

00 

0-54117F 

00 

0,52068E 

00 

AI 1K) = 

0,59386E 

00 

0.6 246 36 

00 

0,633896 

00 

0.62903E 

00 

0,61557E 

00 
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Figure 2. - Blade channel - mean section. 
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Figure 3. - Blade channel - mean section. (Showing typical streamline orthogonal.) 
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Meridional plane wall curvature. 


cttdp|v»NPUT DATA SHEET 

The l^st line in each block displays Roiaan numerals for card sequence and arable numerals for card columns 
vitnin each card. 
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specification of E9.5) usually as £xxx3^zz representing l.xxxaoc multiplied by lO^zz,] 
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Figure 5. - Input data sheet. 







